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Background: The major challenge for nuclear theory is to describe and predict global properties and collective 
modes of atomic nuclei. Of particular interest is the response of the nucleus to a time-dependent external field 
that impacts the low-energy multipole and beta-decay strength. 

Purpose: We propose a method to compute low-lying collective modes in deformed nuclei within the finite 
amplitude method (FAM) based on the quasiparticle random- phase approximation (QRPA). By using the analytic 
property of the response function, we find the QRPA amplitudes by computing the residua of the FAM amplitudes 
by means of a contour integration around the QRPA poles in a complex frequency plane. 

Methods: We use the superfluid nuclear density functional theory with Skyrme energy density functionals, 
FAM-QRPA approach, and the conventional matrix formulation of the QRPA (MQRPA). 

Results: We demonstrate that the complex-energy FAM-QRPA method reproduces low-lying collective states 
obtained within the conventional matrix formulation of the QRPA theory. Illustrative calculations are performed 
for the isoscalar monopole strength in deformed '^^Mg and for low-lying K = quadrupole vibrational modes of 
deformed Yb and Er isotopes. 

Conclusions: The proposed FAM-QRPA approach allows one to efficiently calculate low-lying collective modes 
in spherical and deformed nuclei throughout the entire nuclear landscape, including shape- vibrational excitations, 
pairing vibrational modes, and beta-decay rates. 

PACS numbers: 21.10.Re, 21.60.Jz, 23.20.Js 
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I. INTRODUCTION 

Vibrational modes of atomic nuclei provide crucial in- 
formation about nuclear structure. In particular, collec- 
tive low-lying states contain information about the nu- 
cleonic shell structure, pairing correlations, and nuclear 
deformations [l|, |2[ . Giant resonances tell us about global 
properties of nuclear matter, such as compressibility and 
symmetry energy [SJ, ^ . Electromagnetic strength plays 
an important role in nuclear reactions involving photo- 
nuclear processes, including astrophysical reactions [3, Q 
and the transmutation of nuclear vifaste [7[ . 

The random-phase approximation (RPA) and its su- 
perfluid extension, the quasiparticle random-phase ap- 
proximation (QRPA), are well-established microscopic 
theories describing excitations of many-body systems 
[l|,[2| . QRPA can be viewed as a small-amplitude approx- 
imation of the time-dependent density functional the- 
ory [SJ, |9[. By using nuclear energy density functionals 
(EDFs) applicable to a large portion of the nuclear land- 
scape, a variety of excited modes can be described by 
using QRPA. 

Recently, there has been a considerable progress in the 
area of fully self-consistent QRPA calculations based on 
the nuclear density functional theory. Due to advances 
in high performance computing, deformed QRPA frame- 
works have been developed that can address properties 



of well-bound and loosely-bound nuclei [Tol - [T7| . 

The traditional QRPA methodology is based on a 
generalized eigenvalue problem involving a QRPA ma- 
trix containing the residual two-quasiparticle interaction. 
Because of a large number of quasiparticle states in- 
volved, the dimension of the QRPA matrix is typically 
quite large, especially when spherical symmetry is bro- 
ken. This makes the problem computationally challeng- 
ing. Therefore, in order to reduce the dimension of the 
two-quasiparticle basis, additional cutoffs are imposed on 
the configuration space of QRPA. Such truncations result 
in inconsistencies between the model spaces of Hartree- 
Fock-Bogoliubov (HFB) and QRPA calculations, and can 
result in breaking self-consistency and appearance of spu- 
rious modes [l^ ■ 

To circumvent these problems, efficient methods to 
solve RPA have been formulated in the framework of 
the linear response theory and time-dependent HFB. One 
of these methods is the finite amplitude method (FAM) 
proposed in Ref. [1^. Within FAM, the strength func- 
tion of an arbitrary one-body transition operator can be 
calculated without actually constructing and diagonaliz- 
ing the full (Q)RPA matrix. Instead, the fields induced 
by the one-body transition (driving) operator are cal- 
culated and the linear response problem is solved itera- 
tively. The practical implementation of the FAM requires 
minor extensions to the existing HFB codes to calculate 



the induced fields and, therefore, is fairly straightfor- 
ward. Systematic calculations with the FAM have been 
performed for the electric giant dipole resonances and 
low-lying dipole strength, illustrating computational ad- 
vantages of the method [20, UM ■ The FAM has also been 
extended to the superfluid systems, both spherical [22| 
and deformed [23| . 

Since the FAM equations are solved by introducing a 
small width, an imaginary part of the QRPA frequency, 
the method is very effective for describing excited modes 
in a region of high density of states. However, until now, 
a direct application of FAM to discrete low-lying exci- 
tations has not been fully accomplished. Quite recently, 
an efficient method to evaluate the QRPA matrix using 
FAM has been reported [2J| that significantly reduces 
the computational effort, also enabling computations of 
low-lying discrete QRPA modes. A disadvantage of this 
approach is that a large memory is required to store the 
huge QRPA matrix, which subsequently needs to be di- 
agonalized. 

An alternative technique to solve the linear response 
problem is based on the iterative Arnoldi diagonaliza- 
tion method [23. This method was first implemented 
for spherical systems without pairing and then further 
extended to spherical superfluid nuclei [2g. Because 
the Arnoldi diagonalization algorithm solves the QRPA 
eigenvalue problem in a smaller Krylov-space, the dis- 
crete excitations are within the scope of this method [231 ■ 

The goal of this work is to derive a method to calcu- 
late the discrete low-lying QRPA modes within the FAM 
framework. We shall refer to this new technique as FAM- 
QRPA in the following. Starting from the linear response 
theory, we show in Secs.HIIand lllll that a contour integra- 
tion in the complex frequency plane around a QRPA root 
provides the QRPA eigenvectors. A similar technique 
was proposed in Ref. [23I to solve generalized eigenvalue 
problems. We devise several techniques to compute and 
assess the accuracy of QRPA modes. Next, in Sec. lIVi we 
numerically demonstrate that the discrete FAM- QRPA 
solution for the low-lying states reproduces the modes 
obtained within the conventional matrix formulation of 
QRPA (MQRPA) and we apply FAM-QRPA to collec- 
tive modes in deformed Er and Yb nuclei. Finally, the 
conclusions of our work are given in Sec. |V| 



II. FINITE AMPLITUDE METHOD 

In this section we recapitulate the derivation of the 
FAM equations for superfluid systems following Sec. II 
of Ref. [23. In the FAM formalism, the polarization of 
the system is induced by an external time-dependent field 
F{t) with a frequency w: 



where 



and A = aj^al and B^^ = o.'Lo-l' are two-quasiparticle 
operators. The parameter 77 is a (small) real number to 
expand particle and pair HFB densities to the first order. 
Contrary to Ref. [22| , we assume that F is w-independent 
in all applications in this work. However, our scheme can 
be easily extended to the case where F depends on uj. 

The time-evolution of quasiparticle operators under 
the external field F{t) is determined by the time- 
dependent HFB (TDHFB) equation: 

i^^a^it)^[H{t)+Fit),a^it)], (3) 

where time-dependent oscillation of quasiparticle opera- 
tors is: 



-B„t 



a^(t) = {a^ + (5a^(i)}e*-^" 



(4a) 



6a,it) = 7?^ at {X.^(w)e--* + r,;(c.)e-*} , (4b) 

V 

where i?^ is the one-quasiparticle energy and X^j^^iyj) and 
y^i,(a;) are the FAM amplitudes. 

In terms of time-dependent quasiparticles, the TDHFB 
Hamiltonian can be written as il(t) = Ho+5H{t), where: 



Hq — 2_^ Ei^,B^ 



(5) 



is the HFB Hamiltonian and 

5H{t) = r/ |(5H(a;)e-*'^* + (5Fl"(a;)e'"*} (6) 

with 

5H{u:) ^\Y. {'5^'°(^)^I. + <5<(^)A^.} (7) 

represents a small-amplitude oscillation. 

Inserting ©, ((4|), and © into ^ results in the FAM 
equations: 






_p02 



(8aj 
(8b) 



Ht)^v{ 



^viFe 



■Ft, 



J ' 



(1) 



By expanding 5H'^'^{uj) and (5iJ"^(a;) in terms of X{u!) 
and Y{uj), one obtains: 

l-i' <v' 

+ Y. Bf,^^^,'^>Y^'^'{uj), (9a) 

tJ.'<i^' 



where A and B arc the usual QRPA matrices [3- The 
advantage of the FAM formulation is that the A and 
B matrices do not have to be computed explicitly. By 
substituting © into the FAM equations (|5]), the linear 
response equation becomes: 



A B 



B* A*l''^\0 



1 









(10) 



where TV is the dimension of the two-quasiparticlc space, 
and the matrix elements AT* and F* of X are the QRPA 
amplitudes of the i-th mode with a positive eigenfre- 
quency $7^. There also exists a counterpart QRPA solu- 
tion {Y",X") with a negative eigenfrequency —fti- We 
assume here that all the QRPA frequencies are real, that 
is, S is positive definite. In terms of matrices (|13p . the 
QRPA equation can be expressed as: 



where the sum over two quasiparticle space is restricted 
to quasiparticle indices ^ < v. The FAM equations are 
thus equivalent to the linear response formalism. Fur- 
thermore, the left hand side of PH)) yields the QRPA 
equations when the right-hand side is set to zero. The 
FAM equations ([5]) are solved by using complex frequen- 
cies w-y = w-|-i7, where the imaginary part 7 corresponds 
to a smearing width. 

In terms of the FAM amplitudes X{w^) and Y{uj-y), the 
strength function dB{uj; F)/dio for the operator F can be 
written as: 



dB{uj;F) 1 

= Imb(F;ujy), 



(11) 



5(F; w^) = Y. {^'°*^M^K) + J^^rVK)} ■ (12) 



SX = AfXO, 



(14) 



where O is a 2N x 2N diagonal matrix containing the 
QRPA eigenfrequencies (fii, . . . , JItv, — ^1, • • • , — ^w)- 
The orthonormalization condition for the QRPA eigen- 
vectors is: 



XJ\fX'< = A/". 



(15) 



By applying Eqs. ([Ti|) and ([TSl) . the matrix on the left- 
hand side o f ([T0| can be inverted, which yields the FAM 
amplitudes [2|: 



III. FAM FOR DISCRETE QRPA MODES 



The objective of this work is to formulate a FAM ca- 
pable of describing low-lying discrete QRPA modes. We 
start by introducing the 2N x 2N matrices [2]: 



5 = 



A B 
B* A* 



■^-'J-ll'^ = 



X y 

Y X' 



Yiu;^) 



-i?(a 



r) \p02j 



-x{o - oj^iy^jvx^ 



/^20 
1^02 



(16) 



(13) where I is a 2N x 2N unit matrix and R{uj^) is the 
response function. The explicit form of R{uj~f) is: 



Rfj.i^fi'iy'i^-y) — y ^ 






n^ 


— UJ^ 


Y\ 


X", , 


/if 


/j'l/' 



^,. 









vi^ vi 



(17) 



Substitution of Eq. ()17p into Eq. ([T6l) provides the rela- 
tion between the FAM amplitudes and QRPA amplitudes 




F\^) 



YU^F\0) , X"{0\F\t) 






n. 



■ UJ-y 



n. 



- LO^ 



(18a) 



, (18b) 



where 



(i|i^|0) = ($o|[0.,^]|$o) 

= Y.^KX' + mKl)^ (19) 

(0|F|*) = (<i>o|[()l,F]|<i>o) 

= Y.{f'^:Ixi, + fIIy;,), (20) 

are the QRPA transition strengths between the QRPA 
ground state |0) and i-th excited state |i), |$o} is the 



HFB state, the operator: 



^i ^ ^^ V^/ii/^/ii/ ^/ii/^A'i'J 



(21) 



HKv 



is the QRPA phoiioii operator, and X' and 1?^ are the 
QRPA amplitudes of a state i. 

Equation p^ shows that the FAM amphtudes X{uj^) 
and Y{uj'y) have first-order poles on the real axis at 
io.y = rii and —fti- By calculating the standard FAM 
strength function, approximate positions of the poles of 
the low-lying states of interest can be located. This al- 
lows one to define a closed contour Ci in the complex 
energy plane that encloses the i-th positive pole fi^. Ac- 
cording to Cauchy's integral formula, the contour inte- 
gration of the FAM amplitudes (fT8|) along d gives the 
residue at the i-th pole: 



1 

27ri 



X^,{Lj^)dLu^ = Res{X^,,n,) = XlMF\0), 



(22a) 



27ri 



-- i Y^,{Lu^)du:^ = Res(V,r!,) = y;,(i|F|0). 



(22b) 



The absolute value of the transition strength for the i-th 
QRPA mode can then be expressed as: 



mQ)\^ 



E 
I 1 



1 

27ri 



Xi^,^{ijj^)duj^ 



c. 



Yf_^^{uj^)duj^ 



Ci 



(23) 



where we have used the normalization condition (jisp for 
the QRPA amplitudes. The individual QRPA amplitudes 



AT* and y* can thus be calculated as: 






1 



i(*i^|o)r — f ^M-K)rf^7, 



Y. 



flL' 



2m 



(24a) 
(24b) 



The common phase e'^ = (i|i^|0)/|(i|F|0)| cannot be de- 
termined and remains arbitrary. 

The information about the exact value of the QRPA 
eigenfrequency is not necessary to perform the contour 
integration as long as the corresponding pole is located 
inside the contour. However, it can be calculated from 
the integration of the induced fields. Indeed, from 
Eqs. © and ([13]), one obtains: 



1 



27ri 



— 6 5H^;^{u:,)du:. 



Ci 



2^^ ^ 5Hf,{u:^)duj^ 



(25a) 

{i\F\Q)Y;,{-^r-{Et. + E,)}. 

(25b) 



These 27V equations can be used to compute fi^ , but this 
method is prone to large numerical errors when ampli- 
tudes X* or y* are very small. To this end, a bet- 
ter way of determining the QRPA eigenfrcquencies is 
through an expression derived from Eq. (j25p : 



o2 = ^(|f7,x;,p-|f7,i^, 

I 1 

1 



1 



|(z|F|0)P 



E 



C, 



{{E,, + E„)X^,{u:^) + SHlliij^)) duj, 



2TTi 



C^ 



{{E^ + E,)Y^,{lu^) + SHf^{u^)) du. 



(26) 



The formalism presented above allows one to estab- 
lish an explicit connection between the FAM strength 
function and the smeared QRPA strength function. By 
substituting Eq. p^ into Eq. ([T^ wc obtain: 

' mm? 



S{F,uj,)^-J2 



|(^|F|0)p 



Hi — uj ^ ij Qi + cj + 27 



(27) 



J D 1 



duj 



E 



mm\' 



mm' 



(28) 



According to Eq. ([27]), the discrete QRPA transi- 
tion strength can be directly computed from the FAM 
strength function ([T2]) : 



1 



1(^1^10)1' = ^^ SiF,u)dLo 



(29) 



In summary, as discussed above, there exist several 
techniques, based on the residue at the QRPA pole, to 
calculate discrete transition strengths within the FAM- 
QRPA formalism: 

A: The contour integration of the FAM amplitudes 
X^i,{uj) and Y^^{uj) as in Eq. ([^ : 

B: The contour integration of the FAM strength function 
as in Eq. ([^ : 

C: Individual QRPA amplitudes X' and Y?j^ can be 
found using (j24p to obtain the transition matrix el- 
ement ([T9l): 



D: The QRPA amplitudes X* and F* found with tech- 
nique C are independent of the choice of the external 
field used in FAM-QRPA. Therefore, for example, the 
isoscalar strength associated with the field F' can be 
computed using the QRPA amplitudes obtained in 
FAM-QRPA with the isovector external field F. 



Although ah of these strategies are formally equivalent, 
the technique B is the easiest to implement in the current 
FAM codes. By virtue of D, once the discrete QRPA 
amplitudes have been found for a given state, they can 
be used to calculate a transition matrix element for any 
transition operator. 

If assigned incorrectly, the integration contour C could 
include secondary unwanted poles. (For example, there 
could be two states: a collective one carrying a strong 
transition strength and a nearby-lying non-collective one 
with a negligible contribution to the total transition 
strength.) Since the FAM amplitudes X^i,(w) and Y^i,{ijj) 
p8)) are sums of the residua, the right hand side of 
Eq. ([22]) contains contributions from all the poles in- 
cluded inside C". The calculated transition strength then 
becomes: 



i?(C';F)=^|(z|/'|0)|^ 



(30) 



iec 



Because of the orthogonality of QRPA amplitudes X* 
and y' , the interference terms between different states 
cancel out. Therefore, if C" encircles two or more poles, 
the transition strengths from all those poles contribute 
to the total strength without the interference term when 
techniques A and B are used. Within C, calculated dis- 
crete amplitudes X^^ and F^^ contain a mixture of all 
states inside the contour. However, when applied to 
Eq. ()19|) . the same transition strength as with techniques 
A and B is obtained. However, in the method D, due to 
the incorrect amplitudes X* and F* , the final strength 



i?D(C';f') 



E,6C'(*l^'|o)(*l^|o)' 



E 



iec 



(z|F|0)P 



(31) 



would be incorrect. 

As will be demonstrated in Sec. IIVI we have checked 
numerically that when the contour includes multiple 
poles, techniques A-C indeed yield the total summed 
strength while D does not. This apparent deficiency of D 
can be used to our advantage to verify that the selected 
contour C,; includes only one pole. 

One can also find a clue to correct the assignment of 
the contour by calculating the QRPA eigcnfrequency Vi^ 
using Eq. (^5]) : 

^h - E {w\m?^i) I E K*i^io)i'- (32) 

ieC" iec 

When the contour encloses one collective and one non- 
collective QRPA root with respect to an external field F, 
Eq. (P5| yields the approximate energy of the collective 
state. 



IV. NUMERICAL RESULTS 

To validate the FAM-QRPA formalism discussed in 
the previous section, we carried out numerical compu- 
tations using the FAM framework developed in Ref. [23 



to evaluate the transition strength and the correspond- 
ing residua. Our FAM-QRPA method is based on the 
HFB code HFBTHO [H, H^I, which solves the Skyrme- 
HFB equations in the (transformed) harmonic oscillator 
basis assuming axial and mirror symmetries. The FAM 
equations are solved iteratively by using the modified 
Broyden's procedure (Sll, l32|, which offers a rapid and 
stable convergence, which weakly depends on the magni- 
tude of the imaginary frequency 7. 



A. Test case: monopole strength in ^*Mg 

To compare with full MQRPA, we consider the same 
case of the monopole strength in ^^Mg as discussed in 
Rcf. [ill. We use SLy4 Skyrme EDF [Hi and a con- 
tact volume pairing with a 60 MeV quasiparticle energy 
cutoff and the pairing strength Vb = — 125.20McVfm^ 
for both neutrons and protons. In order to perform exact 
comparison without any truncation at the MQRPA level, 
we take the single-particle basis consisting of N^\^ ~ 5 os- 
cillator shells [23. 

The oblate-deformed HFB minimum of ^^Mg was ob- 
tained at the quadrupole mass deformation /3 = —0.163. 
In this configuration, both neutrons and protons are in 
the superfiuid phase, with pairing gaps An ~ 0.666 McV 
and Ap = 0.654 MeV, respectively. In the FAM calcula- 
tion, we used the value of the parameter r\ = 10~^, which 
was found to provide the best accuracy [23 . For the con- 
vergence criterion of FAM iterations, defined in terms of 
the maximum difference between collective FAM ampli- 
tudes in two consecutive iterations, we used the value of 
10^^; this accuracy is typically reached after about 40 
iterations. 

As for F ^ we consider the isoscalar monopole (ISM) 
and isovector monopole (IVM) operators: 



i=l 

N 



F 



IVM 



A 
eZ 



A Z^'» A 2^' 



i=l 



(33a) 
(33b) 



i=l 



For the integration contours we take circles with radii 
0.02 MeV, centered close to MQRPA frequencies. The 
contour integration is discretized with 11 points, unless 
stated otherwise. 

Figure [1] shows the isoscalar monopole strength func- 
tion at the oblate configuration of ^^Mg calculated with 
the conventional FAM by using three values of 7. The 
strength function obtained with 7 = 0.5 MeV shows a 
very smooth distribution with the broad bumps carrying 
the largest strength. By going to smaller values of 7, one 
reveals the detailed structure of QRPA modes. For exam- 
ple, to separate the smaller first peak at fli = 1.32 MeV 
from the second one at ^2 = 1.37 MeV, a very small 7 - 
of the order of 1 keV - is required. 



TABLE I. Low-lying K = QRPA energies Qt and isoscalar monopole strength |(i|F|0)|^ calculated with MQRPA and FAM- 
QRPA for the oblate configuration of ^*Mg. All the modes with Q,i < 7.5 MeV are listed. The transition strength was computed 
using the techniques A-D described in Sec. IIIII The isoscalar monopole strength FAM-D is calculated from the FAM-QRPA 
amplitudes generated by the external isovector monopole field. The numbers in parentheses denote powers of 10. 



Qi (MeV) 






|(i|^ISM|Q^|2 (g2fjj^4) 






MQRPA FAM 


MQRPA 


FAM-A 


FAM-B 


FAM-C 


FAM-D 


1.3185 1.3183 


5.729(-4) 


5.771(-4) 


5.773(-4) 


5.776(-4) 


5.781(-4) 


1.3731 1.3731 


1.539(-2) 


1.511(-2) 


1.511(-2) 


1.510(-2) 


1.511(-2) 


2.4582 2.4581 


0.1796 


0.1780 


0.1782 


0.1784 


0.1783 


2.5998 2.5975 


2.957(-3) 


3.056(-3) 


3.058(-3) 


3.060(-3) 


3.057(-3) 


3.6687 3.6657 


0.5776 


0.5755 


0.5771 


0.5788 


0.5788 


5.1185 5.1212 


3.539(-4) 


3.744(-4) 


4.040(-4) 


4.360(-4) 


4.345(-4) 


7.4108 7.4084 


0.4900 


0.4820 


0.4834 


0.4848 


0.4848 



TABLE II. Similar as in Table [J but for the isovector monopole modes. 



fli 


(MeV) 






|(i|^IVM|Q^|2 (g2f^4) 






MQRPA 


FAM 


MQRPA 


FAM-A 


FAM-B 


FAM-C 


FAM-D 


1.3185 


1.3183 


1.557(-3) 


1.547(-3) 


1.547(-3) 


1.547(-3) 


1.547(-3) 


1.3731 


1.3731 


5.771(-5) 


5.810(-5) 


5.818(-5) 


5.827(-5) 


5.824(-5) 


2.4582 


2.4581 


1.968(-6) 


1.643(-6) 


1.896(-6) 


2.188(-6) 


2.047(-6) 


2.5998 


2.5975 


8.978(-5) 


8.870(-5) 


8.894(-5) 


8.919(-5) 


8.907(-5) 


3.6687 


3.6657 


1.555(-5) 


8.681(-6) 


1.140(-5) 


1.498(-5) 


1.515(-5) 


5.1185 


5.1212 


3.907(-2) 


3.885(-2) 


3.899(-2) 


3.914(-2) 


3.914(-2) 


7.4108 


7.4084 


1.388(-5) 


2.926(-5) 


2.228(-5) 


1.697(-5) 


1.622(-5) 



> 10^tx10' 



^^^«"-"~*^*^H« 
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FIG. 1. (Color online) The low-lying isoscalar monopole 
strength at the oblate HFB minimum of ^^Mg calculated with 
the conventional FAM-QRPA using three values of smearing 
width 7 (in MeV). 



In Tables U and |TI] wc show the energies and transition 
strength of the low-lying K ~ Q QRPA modes, calcu- 
lated with MQRPA and FAM. Although the centers of 
the contours in the complex w.y-plane are only approxi- 
mately chosen, the low-lying QRPA eigenfrequencies cal- 
culated with FAM are in good agreement with MQRPA 
results. The agreement is excellent for the lowest-lying 
states, which are usually of more interest. The FAM 



transition strength was obtained using the techniques A- 
D described in Sec. IIIII It is gratifying to see that the 
four methods generally agree at least up to two deci- 
mal places, except for the modes carrying very small 
strength (~ 10~'*e^fm'*). The nice agreement between 
FAM-C and FAM-D results demonstrates the consistency 
between the two sets of QRPA amplitudes calculated 
from the isoscalar and isovector external monopole fields. 
The difference between the strengths obtained by the 
MQRPA and FAM is consistent with the convergence 
criteria used in the FAM iterations. 

Figure [2] demonstrates the convergence of the QRPA 
amplitudes against the number of discretization points 
A^disc used in the contour integration. Specifically, it 
shows the orthogonality of the QRPA amplitudes for 
the three lowest QRPA states. The orthogonality be- 
tween the first and third state, and between the second 
and third state, is achieved already at A^disc = 4, while 
the convergence for the pair of first and second states is 
slower. This is not surprising as the energies of the first 
and second QRPA roots differs only by 0.05 MeV; hence, 
and a finer integration mesh is required to remove the 
contribution from the unwanted pole outside of the con- 
tour. Our results show that to obtain the convergence 
for the contour integration, consistent with the accuracy 
required during the regular FAM iterations, taking 11 
points is fully sufficient, at least for the two lowest modes. 

Lastly, wc discuss an example of an incorrect contour 
assignment following the discussion in Sec. IIIII A con- 
tour of radius of 0.2 MeV, centered at 1.3 MeV, includes 
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evaluated with the three-point expression. They are: 
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FIG. 2. (Color online) Convergence of the ortiiogonality of 
tiie states as a number of discretization points of a circular 
contour witii radius 0.02 MeV. Tiie first tiiree low-lying states 
labeled as 1 (1.32 MeV), 2 (1.37 MeV), and 3 (2.46 MeV) are 
shown. 



the first two QRPA solutions at 1.32 MeV and 1.37MeV. 
Such an incorrect choice would be made if the contour 
were determined from the isoscalar strength function cal- 
culated with a resolution of 7 = 0.02 MeV shown in 
Fig. [TJ The calculated isoscalar (isovector) monopole 
transition strength determined according to A-C is 1.57 x 
10~^ (1.61 X 10"'^) e^ fm , which is precisely the sum of 
the two QRPA strengths. However, method D completely 
fails, yielding the values of 6 x 10~* e^ fm (isoscalar) and 
3 X 10~^ e^ fm (isovector), that are clearly off from those 
obtained with the procedures A-C. The QRPA frequency 
calculated from Eq. (|26| with the isoscalar (isovector) 
monopole external field is 1.371 MeV (1.320 MeV). Since 
the isoscalar (isovector) monopole strength of the sec- 
ond (first) QRPA state is larger than that of the first 
(second) QRPA state by two orders of magnitude, the 
QRPA frequency calculated using the contour, which en- 
closes both poles, is close to the energy of the collective 
state. Therefore, the consistency between the results ob- 
tained in methods C and D, together with the value of 
weighted frequency, can be used to find the contour that 
encloses a single QRPA pole. 



Low-lying QRPA modes in deformed rare-earth 
nuclei 



To demonstrate the feasibility of the FAM-QRPA for- 
malism to describe the low-lying collective modes of 
deformed nuclei, we have performed FAM calculations 
for the low-lying K = strength of i66,i68,i72Yb^ and 
^^°Er, which were previously studied with MQRPA in 
Refs. [12, UjI- The calculations were carried out us- 



ing SkM* Skyrme EDF [3^ with the volume pairing. 
The pairing strengths have been adjusted to repro- 
duce the odd-even binding energy difFcrence in -'^^^Yb 



K 



-176 MeV fm"^ and ¥„ 



the quasiparticle energy cutoff Ec 



-218 MeV fm"' 
t = 60 MeV 



for 
and 



Vn = -150 MeV fm"^ and Vp = -177.5 MeV fm'^ for 
-Bcut = 200 MeV. To obtain QRPA amplitudes in FAM, 
we applied the isoscalar quadrupole K — external 
field [i^, and the electric reduced matrix elements B{E2) 
for the excitational modes discussed in Refs. p_2, 13] are 
computed using the technique D described in Sec. IIIII 
The transformed harmonic oscillator basis with 20 major 
oscillator shells was employed. To compute residua, we 
used the circular contours with radii 0.1 MeV for ^^^Yb, 
"^Yb, and ""Er, and with radii 0.01 MeV for ^^SYb. 
The locations of the contour centers estimated from, the 
conventional FAM calculations are: 1.40 MeV, 1.75 MeV, 
1.30 MeV, and 1.30 MeV for ^^^Yh, ^^^Yb, "^Yb, and 
^^°Er, respectively. The results were compared with the 
MQRPA calculations of Refs. [ij, [l^ employing a differ- 
ent HFB solver and an additional cutoff associated with 
the occupation probabilities of canonical states. 

TABLE III. FAM-QRPA energies and B{E2) values of the 
low-lying K = states in ^^^Yb, ^^**Yb, ^^^Yb, and ^''°Er at 
Scut = 200 MeV compared to the MQRPA results of Ref. ^. 
The additional result for '^^^Yb corresponding i^cut = 60 MeV 
is compared to the MQRPA values obtained in Ref. [13 |. 



nucleus 



Q., (MeV) 
MQRPA FAM 



B{E2) (e^b^) 
MQRPA FAM 



166Yb 


1.802 


1.422 


0.0398 


0.0327 


168 Yb 


2.039 


1.747 


0.0343 


0.0186 


172Yb 


1.605 


1.306 


0.0049 


0.0088 


170E, 


1.596 


1.322 


0.0030 


0.0047 


172Yb- 


1.390 


1.319 


0.0050 


0.0092 



a TT" 



60 MeV 



Table IIIII displays the results for excitation energies 
and B{E2) rates of ihc K = Q QRPA modes. For ^^^Yb 
the calculations were carried out with two quasiparti- 
cle cutoffs: -Bcut = 60 MeV and 200 MeV. The FAM- 
QRPA excitation energy is close to the MQRPA value 
with i^cut — 60 MeV [l^, but this agreement does not 
hold when E'cut is increased. Indeed, while our FAM- 
QRPA values weakly depend on -Ecut, a 15% increase of 
the MQRPA energy for "^Yb was reported in Ref. [H 
when going to ^cut = 200 MeV. Interestingly, the B{E2) 
values weakly depend on energy cutoff in both methods. 
However, the B{E2) values obtained in FAM-QRPA are 
twice as large as the MQRPA results. 

For ^^^Yb, ^^^Yb, and ^^"Er, the excitation energies 
obtained in MQRPA are larger by 0.3-0.4 MeV than those 
in FAM-QRPA. The agreement between B{E2) values is 
good in ^^^Yb, but gets worse in the other cases stud- 
ied. It is difficult to speculate what is the origin of those 
differences. We note, however, that (i) the HFB solvers 
used in both calculations are different (see benchmarking 
results in Ref. [35|), and (ii) there are additional canoni- 
cal energy cutoffs in in MQRPA [1^ that are not present 



in our FAM-QRPA method. 



Isoscalar and isovector quadrupole strength (in 
e- fm") of the lowdying K ^0 states in I'^^Yb, ^''^Yb, ^^ 

170t 



TABLE IV. 

„2 f^4\ „f <-V,„ 1^„, !,,;„„ l^ — n r.+ „+^c ;„ 166vk ISSvK i^^Yj-i 

and ^™Er shown in Table [ITI] with Scut = 200 MeV. The 
isoscalar (isovector) strength in FAM-C is calculated using the 
isoscalar (isovector) quadrupole external field. The isoscalar 
(isovector) quadrupole strength in FAM-D is calculated using 
the QRPA amplitudes obtained from the FAM calculation us- 
ing the isovector (isoscalar) quadrupole external field. 



nucleus 



ISQ 

FAM-C FAM-D 



IVQ 
FAM-C FAM-D 



=Yb 
'Yb 
^Yb 
OEr 



299.854 


299.856 


0.585519 


0.585520 


160.126 


160.127 


0.969114 


0.969124 


93.2710 


93.2735 


0.081406 


0.081404 


56.2932 


56.2913 


0.460285 


0.460254 



The isoscalar quadrupole (ISQ) and isovector 
quadrupole (IVQ) strengths are displayed in Table IIVI 
for the deformed nuclei shown in Table Hill The strengths 
obtained from the QRPA amplitudes derived from two 
external quadrupole fields agree excellently. This result 
clearly shows that the contours used in these calculations 
enclose only a single QRPA pole. 



formulation is more appropriate. 

The FAM-QRPA method has been benchmarked and 
tested by comparing with MQRPA results for an oblate 
configuration of ^^Mg. Illustrative examples of large- 
scale calculations have been presented for the K = 
isoscalar and isovector quadrupole modes of deformed 
rare-earth nuclei i^^^Yb, i^Syb, i72Yb, and i™Er. 

Our results demonstrate that the proposed formula- 
tion of FAM-QRPA can be used as an efficient tool to 
calculate discrete QRPA modes of heavy, deformed, and 
superfluid nuclei. Once the contour around the mode of 
interest is specified, the FAM-QRPA method allows one 
to perform a fully self-consistent QRPA calculation em- 
ploying the same model space as in HFB. Thanks to the 
rapid convergence achieved with Broyden's method used 
in our implementation, FAM-QRPA is amenable to high- 
performance parallel computing. This offers promise of 
systematic calculations of various kinds of low-lying exci- 
tations and decays over the entire nuclear landscape. Of 
particular importance are QRPA studies of low-energy 
dipole and quadrupole states, (3 decays, and (3/3 decays. 
The work on extending the FAM-QRPA formalism to 
K ^ and charge- exchange modes is in progress. 
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